%% add paths
%% Non-selective 
clear all
addpath(genpath('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\'))

clear PD_NSRF4; clear FlipAngle_NSRF4;
% Slice 1: calculating mean signal intensity for different vials 
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip5.mat')
PD_NSRF4 = abs(TimeAverageCombined{1});
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip10.mat')
FlipAngle_NSRF4(:,:,1) = abs(TimeAverageCombined{1})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip15.mat')
FlipAngle_NSRF4(:,:,2) = abs(TimeAverageCombined{1})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip20.mat')
FlipAngle_NSRF4(:,:,3) = abs(TimeAverageCombined{1})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip25.mat')
FlipAngle_NSRF4(:,:,4) = abs(TimeAverageCombined{1})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip30.mat')
FlipAngle_NSRF4(:,:,5) = abs(TimeAverageCombined{1})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip35.mat')
FlipAngle_NSRF4(:,:,6) = abs(TimeAverageCombined{1})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip40.mat')
FlipAngle_NSRF4(:,:,7) = abs(TimeAverageCombined{1})./PD_NSRF4;

load('mask_slc1.mat')
mask = mask(:,:,2:10);
clear Signal_NSRF4_slc1
for jj =1:size(FlipAngle_NSRF4,3)
    for mm = 1:size(mask,3)
         temp1 = circshift(imrotate(abs(FlipAngle_NSRF4),-90),[0 1]);
         temp = (mask(:,:,mm).*temp1(:,:,jj));
         Signal_NSRF4_slc1(mm,jj)= sum(temp(:))./sum(sum(mask(:,:,mm)));
    end
end

clear PD_NSRF4; clear FlipAngle_NSRF4;
% Slice 2: calculating mean signal intensity for different vials 
addpath(genpath('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\'))
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip5.mat')
PD_NSRF4 = abs(TimeAverageCombined{2});
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip10.mat')
FlipAngle_NSRF4(:,:,1) = abs(TimeAverageCombined{2})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip15.mat')
FlipAngle_NSRF4(:,:,2) = abs(TimeAverageCombined{2})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip20.mat')
FlipAngle_NSRF4(:,:,3) = abs(TimeAverageCombined{2})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip25.mat')
FlipAngle_NSRF4(:,:,4) = abs(TimeAverageCombined{2})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip30.mat')
FlipAngle_NSRF4(:,:,5) = abs(TimeAverageCombined{2})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip35.mat')
FlipAngle_NSRF4(:,:,6) = abs(TimeAverageCombined{2})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip40.mat')
FlipAngle_NSRF4(:,:,7) = abs(TimeAverageCombined{2})./PD_NSRF4;

load('mask_slc2.mat')
mask = mask(:,:,2:10);
clear Signal_NSRF4_slc2
for jj =1:size(FlipAngle_NSRF4,3)
    for mm = 1:size(mask,3)
         temp1 = circshift(imrotate(abs(FlipAngle_NSRF4),-90),[0 1]);
         temp = (mask(:,:,mm).*temp1(:,:,jj));
         Signal_NSRF4_slc2(mm,jj)= sum(temp(:))./sum(sum(mask(:,:,mm)));
    end
end


% Slice 3: calculating mean signal intensity for different vials 
clear PD_NSRF4; clear FlipAngle_NSRF4;
addpath(genpath('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\'))
load('mask_slc2.mat')
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip5.mat')
PD_NSRF4 = abs(TimeAverageCombined{3});
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip10.mat')
FlipAngle_NSRF4(:,:,1) = abs(TimeAverageCombined{3})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip15.mat')
FlipAngle_NSRF4(:,:,2) = abs(TimeAverageCombined{3})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip20.mat')
FlipAngle_NSRF4(:,:,3) = abs(TimeAverageCombined{3})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip25.mat')
FlipAngle_NSRF4(:,:,4) = abs(TimeAverageCombined{3})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip30.mat')
FlipAngle_NSRF4(:,:,5) = abs(TimeAverageCombined{3})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip35.mat')
FlipAngle_NSRF4(:,:,6) = abs(TimeAverageCombined{3})./PD_NSRF4;
load('H:\SVN\hassan-papers\WMAMP_Perfusion\CodeFig3\Data\Non_Selective_Flip40.mat')
FlipAngle_NSRF4(:,:,7) = abs(TimeAverageCombined{3})./PD_NSRF4;

load('mask_slc3.mat')
mask = mask(:,:,2:10);
clear Signal_NSRF4_slc3
for jj =1:size(FlipAngle_NSRF4,3)
    for mm = 1:size(mask,3)
         temp1 = circshift(imrotate(abs(FlipAngle_NSRF4),-90),[0 1]);
         temp = (mask(:,:,mm).*temp1(:,:,jj));
         Signal_NSRF4_slc3(mm,jj)= sum(temp(:))./sum(sum(mask(:,:,mm)));
    end
end


clear T1_phantom
T1_phantom(1) = 354;
T1_phantom(2) = 415;
T1_phantom(3) = 498;
T1_phantom(4) = 515;
T1_phantom(5) = 621;
T1_phantom(6) = 724;
T1_phantom(7) = 922;
T1_phantom(8) = 943;
T1_phantom(9) = 1178;


clear alpha
alpha(1) = 10;
alpha(2) = 15;
alpha(3) = 20;
alpha(4) = 25;
alpha(5) = 30;
alpha(6) = 35;
alpha(7) = 40;

ColorUse{1} = '-Ok'; 
ColorUse{2} = '-Or';
ColorUse{3} = '-Ob';
ColorUse{4} = '-Og';
ColorUse{5} = '-Om';
ColorUse{6} = '--Ok';
ColorUse{7} = '--Or';
ColorUse{8} = '--Ob';
ColorUse{9} = '--Og';

%% Calculating myocardial contrast 
% myocardial contrast = ration of mean signal intensity for peak myocardial contrast vial to other 8 vials

figure; 
for jj = 2:size(mask,3)
MC_NSRF4_slc1(:,jj-1) = squeeze(abs(Signal_NSRF4_slc1(1,:))./(abs(Signal_NSRF4_slc1(jj,:))));
plot([10:5:40],MC_NSRF4_slc1(:,jj-1),ColorUse{jj-1},'LineWidth',4,'MarkerSize',10);
hold on;
end
axis square; set(gca,'LineWidth',5,'tickLength',[0 0],'FontSize',30,'FontWeight','bold')
ylim([0 2.5])
xlim([10 40])


figure; 
for jj = 2:size(mask,3)
MC_NSRF4_slc2(:,jj-1) = squeeze(abs(Signal_NSRF4_slc2(1,:))./(abs(Signal_NSRF4_slc2(jj,:))));
plot([10:5:40],MC_NSRF4_slc2(:,jj-1),ColorUse{jj-1},'LineWidth',4,'MarkerSize',10);
hold on;
end
axis square; set(gca,'LineWidth',5,'tickLength',[0 0],'FontSize',30,'FontWeight','bold')
ylim([0 2.5])
xlim([10 40])


figure; 
for jj = 2:size(mask,3)
MC_NSRF4_slc3(:,jj-1) = squeeze(abs(Signal_NSRF4_slc3(1,:))./(abs(Signal_NSRF4_slc3(jj,:))));
plot([10:5:40],MC_NSRF4_slc3(:,jj-1),ColorUse{jj-1},'LineWidth',4,'MarkerSize',10);
hold on;
end
axis square; set(gca,'LineWidth',5,'tickLength',[0 0],'FontSize',30,'FontWeight','bold')
ylim([0 2.5])
xlim([10 40])

%% Difference maps for myocardial contrast for different slices
clear str1
figure; 
for jj = 1:size(MC_NSRF4_slc3,2)
plot([10:5:40],abs(MC_NSRF4_slc1(:,jj)-MC_NSRF4_slc2(:,jj)),ColorUse{jj},'LineWidth',4,'MarkerSize',10);
hold on;
end
% legend(str1)
axis square; set(gca,'LineWidth',5,'tickLength',[0 0],'FontSize',30,'FontWeight','bold')
ylim([0 0.1])
xlim([10 40])

clear str1
figure; 
for jj = 1:size(MC_NSRF4_slc3,2)
plot([10:5:40],abs(MC_NSRF4_slc1(:,jj)-MC_NSRF4_slc3(:,jj)),ColorUse{jj},'LineWidth',4,'MarkerSize',10);
hold on;
end
% legend(str1)
axis square; set(gca,'LineWidth',5,'tickLength',[0 0],'FontSize',30,'FontWeight','bold')
ylim([0 0.1])
xlim([10 40])


clear str1
figure; 
for jj = 1:size(MC_NSRF4_slc3,2)
plot([10:5:40],abs(MC_NSRF4_slc2(:,jj)-MC_NSRF4_slc3(:,jj)),ColorUse{jj},'LineWidth',4,'MarkerSize',10);
hold on;
end
% legend(str1)
axis square; set(gca,'LineWidth',5,'tickLength',[0 0],'FontSize',30,'FontWeight','bold')
ylim([0 0.1])
xlim([10 40])
